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This report outlines a computer program for calculating, directly 
from the atmospheric data, the location of ground level acoustic foci 
resulting from meteorological conditions in a three-layered atmosphere. 
Also, a procedure is given whereby one obtains the average intensity 
level in the immediate neighborhood of a focus. Finally, it is shown, 
by means of examples, that a systematic use of the computer program 
can lead to a deepened insight into the relationship of the focus 
location and the atmospheric parameters. 
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FOCAL POINT COMPUTATION IN A THREE-LAYERED ATMOSPHERE 

SUMMARY 


This report outlines a computer program for calculating, directly 
from the atmospheric data, the location of ground level acoustic foci 
resulting from meteorological conditions in a three-layered atmosphere. 
Also, a procedure is given whereby one obtains the average intensity 
level in the immediate neighborhood of a focus. Finally, it is shown, 
by means of examples, that a systematic use of the computer program 
can lead to a deepened insight into the relationship of the focus 
location and the atmospheric parameters* 


I, INTRODUCTION 


As a result of the procedures used in obtaining atmospheric data, 
a velocity-of-sound profile is normally represented by a many-sided 
polygon, where the velocity is customarily plotted as a function of 
height. This polygon, corresponding to a multi-layered atmosphere, is 
assumed to be a more or less accurate approximation of the true vertical 
velocity distribution. From a given profile one obtains velocities and 
heights which, in turn, are used as input to the computer program written 
to calculate sound ray patterns and sound ray landing distances. Normally, 
the only output of concern is the set of landing distances associated 
with a given velocity profile, since, from this set, a determination may 
be made as to the existence (or nonexistence) of one or more focal points. 
This examination procedure is both laborious and time consuming, and 
sometimes fails even to indicate the presence of a focus. 

It would, therefore, be advantageous to have some rapid method not 
only of calculating focal distances directly from the atmospheric data 
but also of insuring that all focal points have been obtained. Such a 
direct method would be very useful in handling several of the problem 
areas which exist in sound propagation studies, resulting from the vast 
amount of output generated by present computer programs. For example, 
consider the problem of determining the effect on focal distances of 
variations in the atmospheric parameters. It is almost impossible, using 
present techniques, to make any sort of comprehensive survey of the 
relationship between the atmospheric parameters and focal distances 
because of the overwhelming volume of computer output which would have 
to be analyzed. Let us consider another problem: Assume that 
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there exists a profile which, contrary to observation, indicates no focus 
formation when analyzed by existing techniques. One might wish to explore 
the possibilities of varying the input parameters in various ways calcu- 
lated to result in a profile which does give focus formation. (However, 
one is entitled to vary the atmospheric parameters only within the accuracy 
limits of the instruments which make the meteorological measurements.) 

This problem can be dealt with, in a reasonable length of time, only by 
an examination of a number of different velocity profiles, each of which 
is only slightly different from the original. In addition, one could 
examine the possibility of varying a profile in some manner in order to 
obtain a focus at a given location. 

A direct method of focal point computation would be of significant 
use in handling these problems. It could also be useful when one is 
faced with a problem which requires analyzing a large amount of computer 
output. However, from both a mathematical and a numerical point of view, 
such a method is nearly impossible to obtain in the case of a general 
multi-layered atmosphere. This is explained as follows. 

The condition for a focal point is that the equation 
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should be satisfied.* (The D's and 5's are functions of the atmospheric 
parameters, k is the number of layers in the atmosphere, and q is related 
to the focal ray's angle of departure, 0 q, which is the unknown parameter.) 
If a q can be found which satisfies equation (1), it can then be used to 
calculate the focal distance, Xp, from 
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Therefore, the problem of direct focal point computation hinges on 
the ability to solve equation (1), for q. 

When k = 2, corresponding to a two-layered atmosphere, equation (1) 
can be solved quite readily (cf. Ref, 1, pp. 30-39). 

When k = 3, corresponding to a three-layered atmosphere, equation (1) 
is not quite so readily solved for q. For a few values of the D's and 
6’s, analytic solutions exist, but, in general, q must be found by some 
numerical methods, the best of which were found to be iterative procedures 
A complete mathematical treatment of equation (1), for k = 3, may be 
found in Reference 2. 

When k = 4, corresponding to a four-layered atmosphere, equation (1) 
can be rewritten as a twelfth-degree polynomial in q. This polynomial 
not only has no analytical solution, but also is not readily amenable 
to numerical root-finding techniques because the coefficients are 
ill-conditioned. When k > 4, expression (1) becomes even more difficult 
to handle. 

The mathematical difficulties associated with the multi-layered 
case prohibit the development of a direct method of focal point calcu- 
lation, except in the two simple cases where k = 2 and k = 3. Fortunately 
many multi-layered situations can be approximated by these simple cases 
when information is sought as to possible focus formation. This is not 
to say that the simple cases should be substituted for the more complex 
cases every time, but, in general, it can be said that if the two- or 
three-layered approximant provides a focus, the multi-layered case will 
also.* 

Therefore, it was decided to take the mathematical techniques 
presented in Reference 2 and to develop a computer program for direct 
focal point computation in the three-layered atmosphere. With the aid 
of such a program, it would be possible to gain a deepened insight 
into the relationship between atmospheric parameters and focus formation. 

This report, then, describes the computer program and how it can 
be used. Also given are some preliminary results of investigations 
carried out in some of the areas discussed previously. A copy of the 
program, written for the CDC 3200 computer, by Mr. J. Hilliard, is 
included as Appendix I for those who wish to explore its possibilities 
further. 


* Reference 5 also points this out. 


II. COMPUTER PROGRAM 


To facilitate discussion of the computer program, a sample page of 
computer printout (Table 1) is included here. In general, the computer 
program allows one to input any set of initial velocities and heights 
(Table 1, Lines 9-10), any set of increments greater than 1 x 10~ 4 to 
be applied to these initial values (Table 1, Lines 12-13), and any set 
of integers, from 1 to 999 (Table 1, Line 7), corresponding to the number 
of times each increment is to be applied. One other quantity, tolerance 
1, (Table 1, Line 4), must be specified as a criterion for the convergence 
of an iteration procedure. The example shows that only the initial 
velocity, V 1? will be incremented, by an amount 0.0125 a total of 35 
times, with a convergence criterion of 1 x 10~ 7 . 

Furthermore, the user has control (by means of two sense switches) 
of what information is printed out: Sense switch 1 gives the option of 
printing each step of an iteration (cf. Table 1, Lines 35-49); sense 
switch 2 gives the option of printing the quantities, exemplified by 
Table 1, Lines 22-25, computed from the velocities and heights. Normally, 
the sense switches are off since, if one is running a particularly large 
number of cases, there is no wish to obtain any more printout than 
necessary. However, either sense switch may be turned on during execution 
of the program if the user desires to monitor the progress of the program, 
or if he desires to verify that everything is being computed properly. 

Although the computer program is basically designed to determine foci 
resulting from atmospheric conditions in the third layer, it also will 
determine whether or not there are foci resulting from atmospheric condi- 
tions in the second layer. This section was included not only for 
completeness but also to take care of the following eventuality. Suppose 
that in the course of a systematic survey, where several velocities and 
heights are being varied simultaneously, it happens that a certain combina- 
tion of the parameters results in two adjacent atmospheric layers having 
the same gradient.* When this occurs, one no longer has a three -layered 
atmosphere but has, instead, an atmosphere of only two layers, which must 
be handled differently. 

Therefore, there has been built into the program existence criteria 
for foci to result from atmospheric conditions in both the two-layered 
and three- layered atmospheres. 


* The gradient is defined as the rate of change of propagation velocity 
of sound with respect to height. 
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FOCAL CONDITION* \ 8 . 7311491370 E - 11 
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XFs 9 . 1184250 7 7 IE 03 SCNDD* 4 . 0770965498 E 05 

OMEGA* 9 . 3076936468 E -09 All* 7 . 7389463162 E -0 1 
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There are two headings under which information is printed concerning 
the results of testing the atmosphere against the foci existence criteria; 
"FOCUS FROM SECOND LAYER" (Table 1, Line 28) and "FOCUS FROM THIRD LAYER" 
(Table 1, Line 33). If the existence conditions are not satisfied in 
either or both layers, the message f, N0 FOCUS" is printed (cf. Table 1, 

Line 30). If, however, the existence conditions are satisfied, meaning 
that a focus does exist, the focal distance, along with other pertinent 
information, is printed. On the occasions when there are two foci 
resulting from the third layer, the information is printed out in two 
consecutive sets. In the example (Table 1) there is no focus from the 
second layer and only one from the third. 

There are nine quantities, all related to the solution of equation (1) 
or to xp, which are printed out each time the program determines that there 
is a focus, whether resulting from atmospheric conditions in the second 
or third layer. These nine quantities are discussed in detail as follows. 

1. FOCAL CONDITION (Table 1, Line 51) - This quantity is a measure 
of how accurately q has been determined, and should be very close to zero, 
since it is an evaluation of equation (1), using the value of q which has 
been obtained. In the example, FOCAL CONDITION ~ 10" 11 , which leads one 
to believe q is quite accurate. 

2. Q (Table 1, Line 52) - This is the solution of equation (1), no 
matter how obtained. 

3. CTST (Table 1, Line 52) - This quantity is cotangent 0^, where 0* 
is the departure angle, measured from the horizontal, of the focal ray, 
when leaving the sound source. 

4. THOS (Table 1, Line 53) - This is 0^, having dimension degrees. 

5. THETM (Table 1, Line 53) - This quantity is 0 O ma ximum> which is 
the largest angle of departure for rays which return to the sound source 
horizontal, y 0 , from the third layer.* It is necessary that THETM be 
computed and printed out because of the mathematical treatment of equation 
(1). All the iteration procedures of Reference 2 are derived under the 
assumption that the third layer is infinitely extended, although, in 
practice, this is not the case. Hence, one must check the value (THOS - 
THETM). If this difference is negative, the focal ray is returned from 
the third layer. A positive difference means that the focal ray traverses 
the third layer entirely and enters the fourth, the third layer being too 
thin in such a case; the presumptive focal ray does not return from it and 


* If the focus results from the second layer atmospheric conditions, THETM 
is the maximum angle of departure for rays which return to y Q from the 
second layer. 
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focus formation is precluded. In the example, (THOS - THETM) < 0; the 
third layer is extended far enough to enable the focal ray to turn around 
and return to the sound source horizontal where it gives rise to a focus. 

6. XF (Table 1, Line 54) - This quantity is the focal distance, 
counted from the sound source. It has dimension meters if [v] = m/sec, 

[y] = meters; or feet if [v] = ft/sec, [y] - ft. 

7. SCNDD (Table 1, Line 54) - This quantity is the second derivative 
of x s with respect to 0 O , evaluated at 0 O = 9* and is printed out as a 
check.* It is used later in the computation of OMEGA and AIL. 

8. OMEGA (Table 1, Line 55) - This is a quantity which is indicative 
of how strongly the focus intensity becomes infinite, and is proportional 
to the intensity average in a differential vicinity of a focus. 

9. AIL (Table 1, Line 55) - This quantity is ( [OMEGA | • Xp 2 ). The 

average intensity in the vicinity of a focal point is found by I = AIL/ e 
where e is an angular difference, taken at the sound source, between the 
focal ray and its differential neighbor. This average intensity can then 
be compared with some reference intensity, for instance, the average 
intensity associated with a homogeneous medium. The choice of the magni- 
tude of e is left to the discretion of the user, the only stipulation 
being that it should be chosen very small (on the order of 10' 3 or less). 
This is because the expression for cj has been derived for £ 0. Also, 

the intensity decreases very rapidly if one moves from a focus; one should 
stay very close to the focus to obtain a meaningful average. Therefore, 

£ should be chosen small enough in order for I to be larger than the 
reference value. If this reference value is taken as the average intensity 
I* near a standard focus, one usually chooses the same value of £ for the 
two foci, so that £ cancels when forming the intensity level 10 log 1Q I/I*. 

With this brief explanation of what the computer generates, the 
remainder of this report will deal with two problems which can be handled 
quite easily. 


III. OBTAINING A FOCUS AT A DESIRED LOCATION 


One of the problems mentioned earlier is that of varying the velocities 
and/or heights of a given velocity polygon in various manners calculated to 
result in focus formation at a particular location. The only way this can 
be handled is to run a number of closely related polygons through a computer 




* Since dx s /d0 Q = 0 at a focal point, SCNDD is identically equal to the 
curvature at xp. 
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program for calculating sound ray landing distances, examine these landing 
distances and see if any of them are the desired focus. If not, more 
variations must be made and more polygons run. The time required to do 
this and the amount of computer output to be examined are dependent on 
how "skillful 11 the person running the program is in making the variations. 
There is no guarantee that any polygon can be found which yields a focus 
at a desired location. 

The author of this report made such a study (cf. Ref. 3) in an 
attempt to obtain a focus at a distance of 27.2 km, by varying the 
parameters in an eleven- layered atmosphere. Observation suggested the 
presence of this focus, but it could not be obtained from the original 
velocity profile. This study took two weeks and an examination of 
fifty-six different sets of sound ray landing distances before one was 
found to yield a focus at 27.2 km. 

Since this time, the program for determining the existence and loca- 
tion of focal points resulting from atmospheric conditions in a three- 
layered medium has been developed. It was decided to make the same study, 
but now replacing the eleven-layered atmosphere by one of only three 
layers and varying the parameters within the measurement accuracies. 

Figure 1, solid line, is a plot of the original eleven-layered 
atmosphere, with the dashed line being its initial three -layer approximant. 
Table 2 gives a summary of the input parameters to the focal point deter- 
mination program. The initial velocities and heights are given in columns 
1 and 2, respectively, column 3 gives the increments of the velocities 
with column 4 denoting the number of times each is to be applied. This 
table shows that only two of the possible eight parameters were varied. 


TABLE II 


VI 

(m/sec) 

YI 

(m) 

AVI 

(m/sec) 

CVI 

339.98 

0 

0.0 

1 

345.54 

721 

0.0 

1 

345.25 

858 

0.0125 

50 

345.89 

1016 

0.0125 

24 
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338 340 342 344 346 348 

V (m/sec) 


FIG. 1. VELOCITY OF SOUND POLYGON 
ACTUAL VERSUS THREE LAYER APPROXIMANT 




These increments resulted in a total of twelve hundred different 
profiles which were tested one by one for possible focus formation. The 
total computer running time was 28 minutes, which is considerably less 
than the more than three hours of running time the original study required, 
even though nearly twenty- two times as many profiles were examined. 
Moreover, the three-plus hours does not include the time required to check 
each profile for the focus, whereas the 28-minutes does. The study in 
Reference 3 did not include any estimate of the average intensity level 
in the vicinity of the focus, whereas the three-layered atmosphere program 
includes a measure of the average intensity level in the neighborhood of 
a focus. (See Appendix II for the derivation of this average intensity 
level . ) 

The results of this new study are summarized in Figures 2 and 3. 

Figure 2 gives the focal distance, Xp, plotted parametrically as a function 
of V 2 and V 3 . Immediately one sees that if he is searching for a focus 
at 27.2 km, he has a whole host of (V 2 , V 3 ) combinations from which to 
choose. Since the computer program calculates an oj for each x F , the 
average intensity level in the vicinity of each xp can be obtained. 

These results are given graphically in Figure 3, where the intensity 
level is plotted parametrically as a function of V 2 and V 3 .* (The 
average intensity levels were computed under the assumption of D = 210 db 
and £ = 10" 4 ; cf. Appendix II.) 

As an example, let V 2 = 345.5 m/sec. Figure 2 shows that xp = 27.2 
km crosses V 2 = 345.5 m/sec at V 3 ~ 346.09 m/sec. For the same V 2 and 
the resulting V 3 , Figure 3 indicates an average intensity level of 122 db. 
This same Figure 2 can be used to determine whether or not a particular 
(V 2 , V 3 ) combination would yield a focus. When a focus does result for 
any (V 25 V 3 ) combination, Figure 3 will indicate the average intensity 
level about this focus. 


IV. VARIATION OF X p RESULTING FROM VARIATIONS IN V G AND V 3 

As another example of the versatility of the focal point computation 
program, the reader's attention is called to Table III. 


* In both Figures 2 and 3, each continuous curve represents a constant 
value of V 3 , with V 2 being the abscissa. 


10 



Atmospheric Conditions 


(m) 


V 0 « 339.98 
V 1 * 345.54 


V 3 {m/sec) *345.8900 


345.9025 


345.9150 


KMP 


' / 345.9400 y 


r / 345.9525 s 




345.2500 .3125 .3750 .4375 .5000 .5625 .6250 .6875 .7500 345.8125 

V 2 (m/sec) 

FIG. 2. VARIATION IN x F RESULTING FROM 
VARIATIONS IN V 2 AND V 3 
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TABLE III 

Input Parameters for Survey 


^initial 

(m/sec) 

^initial 

(m) 

AV 

(m/sec) 

AY 

On) 

CVI 

CYI 

335.00 

0.0 

0.0625 

0.0 

21 

1 

336.00 

100.0 

0.0 

0.0 

1 

1 

336.25 

125.0 

0.0 

0.0 

1 

1 

336.25 

170.0 

0.0625 

0.0 

37 

1 


Tolerance 1 = 1.0 3 

K 10‘ 8 




This table represents the input to the program for the purpose of 
making a survey of what happens to the focal distance, xp, if V Q and V 3 
are allowed to vary simultaneously. These particular values were chosen 
because various combinations result in zero, one or two focal points. 

Seven hundred and seventy-seven different (V 0 , V 3 ) combinations were 
examined, the results being summarized in Table IV.* A glance at this 
table shows immediately which (V Q , V 3 ) combinations yield no focal points, 
which combinations yield one focal point, and, perhaps most significantly, 
which combinations yield two focal points. Table IV also serves to stress 
a point which cannot be overemphasized; that is, that small changes in 
velocities may make significant changes not only in focus location but 
also in their very existence. For illustration, examine V 3 = 336.5 m/sec. 
For 335.0 m/sec ^ V c ^ 336.125 m/sec, there is no focus formation, and 
for 336.25 m/sec ^ V Q ^ 336.4375 m/sec, there is only one focus formed. 

But for V Q = 336.1875 m/sec there are two foci formed by rays returning 
from the third layer. Hence, for a very limited range of V Q , i.e., from 
V Q = 336.125 to V 0 = 336.25, three different focus formation situations 
arise. The conclusion must be that V Q , in this region, has to be known 
quite accurately. 

To further illustrate the amount of information which can be gained 
from proper use of the computer program, two additional figures are 
included here: Figure 4, which is xp plotted parametrically as a function 


A 

* 0 implies no focal point, 1 implies one focal point, 2 implies two 
focal points. 
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TABLE IV 



337.4375 


























































x F (km) 



5.0 -F 1 1 1 1 h- 

335.0 335.5 336.0 336.5 337.0 337.5 


Vq (m/sec) 

FIG. 4. Xp VERSUS V Q , WITH V 3 CONSTANT 
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of V Q and V 3 and Figure 5, which is the average intensity level plotted 
as a function of the same arguments, V c and V 3 .* Only a few of the 37 
different V 3 curves are shown. 

A close examination of Figure 4 reveals several interesting features, 
not the least of which is double-focus formation from the third layer, 
cf. the two -branched curve V 3 ® 336.6875 m/sec. Note, however, that for 
V Q large enough, in this case V Q > ~ 336.1825 m/sec, only one focus is 
formed. Further, for V Q small enough, V 0 < ^ 335.25 m/sec, no focus 
whatsoever is formed. As one would expect from an examination of the 
two-branched curve in Figure 4, there is a value of V Q (« 335.75) for 
which the two focal points converge. This focal point convergence leads 
to inordinately high intensity levels, as the V 3 = 336.6875 m/sec curve 
in Figure 5 indicates. It has been shown analytically that the two 
branches actually converge to an infinite value of the average intensity 
level, so that the focus is uncommonly strong. 

Furthermore, the curve labeled "second layer focal line" in Figure 4 
denotes the focal points obtained as a result of atmospheric conditions 
in the second layer. They approach within 200 meters the cusp focus of 
the two-branched V 3 = 336.6875 m/sec curve. For this particular set of 
atmospheric conditions, the neighborhood around 6880 meters appears as 
one of very high intensity level. An examination of the corresponding 
average intensity level curves of Figure 5 verifies this conclusion. 

(When determining the average intensity level for a particular V Q value, 
the upper branch of xp, when V 3 = 336.6875 m/sec, is associated with the 
lower branch of the average intensity level curve.) 

Another undesirable situation is typified by the V 3 = 336.75 m/sec 
curve in Figure 4. This particular V 3 -value permits only one focal point 
to develop from the third layer. However, the second layer focal line 
crosses this curve at V Q ~ 335.4375 m/sec. Here is a case of identical 
xp's occurring from two different layers, which should lead to very high 
average intensity levels. 


CONCLUSIONS 


The computer program for direct computation of two- and three-layer 
focal points facilitates and expedites, to a significant degree, the 
systematic study of the relationship between atmospheric parameters and 
ground level acoustic focus location. It cannot answer all questions, 
neither does it solve all problems concerning focal points; but it does 
provide a valuable tool for further investigations both in the theoretical 
field and in the area of testing a given meteorological situation for 
possible focus formation. 


* The average intensity level is computed under the assumptions D = 210 db 
and e = 1.0 x 10 “ 4 . 


17 



APPENDIX I 


PROGRAM MABRY 

1000 R f ad ( 60 * 1 )V0l*VlJtV2I»V3I»Y0l»YJ.I»Y2l»Y3I 

RFAn(60*l?00)DVO»DVI»DV2*DV3*DYO*DYl»DY2tDY3 
1200 FORMAT(BFR.O) 

1 FORMAT (4F15.0/4F15.0) 

RFAD(60*6 )N0*N1 ,N2 »N3 *N4 »N5 »N6»N7 

6 FORMAT (815) 

RFAD(60*?)T0L1 *LAST 

2 FORMAT(F10.0*T5) 

WRITE(61»3)T0L1 

3 FORMATC 1H1 »49X,33HTHRFF LAYFRFD ATMOSPHFRF ANALYSTS///54X * 1 3HTOLER 
1ANCE1- F15.8) 

WRITF(61 »40)N0,N1 »N2 *N3 »N4 ,N5»N6 *N7 

40 FORMAT(//12X» 5HCV0= I3»6X»5HCV1= I3*6X»5HCV2= I3,6X»5HCV3= I3,6X»5 
1HCY0= 1 3 * 6X * 5HCY1 = 1.3 ,6X , 5HCY? = I3,6X»5HCY3= 17) 

WRITF(61»A1 1 VOT »V1 1 *V2T * V? T *Y0I *Y1T » Y? T » Y3 1 

41 FORMAT( /10X*5HV0I= F15.10»10X*5HV1 1= FI 3 . 10 * 1 OX , FmV 2 I = F15.10*10X» 
15HV3I= Fl 5.10/10X»5HY0T= Fl 5. 1 0 » 1 OX ♦ 5HY1 I = Fl 5 . 1 0 , 1 OX » 5HY? I = FJ5.1 
20 » 1 OX » 5HY3 1 = FJ5.10) 

WR I TF ( 61 ,42)nV0»DVl»nV2»nV3*r>Y0»nYl *nY2»r>Y3 
4? FORMAT C /10X,5HDV0= F15 . 1 0 » 1 OX , 5HPV1 = Fl 5 . 1 0 » 10X , 5HDV2= F15.10*10X» 
15HDV3= F15.10/10X,5HPY0= F15.10»10X» 3HRY1 = F15. 10 ,10X» 5HDY2= F15.1 
20*1 OX »5HRY3 = F15.10) 

V0=V0I 

DO 30 10=1 »NO 
V1=V1I 

DO 31 11=1, Ml 
V2 = V2 I 

DO 32 12=1 »N2 
V3=V3 I 

DO 33 1 3 = 1 »N3 
Y0=Y0I 

DO 34 14=1, N4 
Y1=Y1I 

DO 35 15=1, N5 
Y2=Y2! 

DO 36 16=1, N6 
Y3= Y3 I 

DO 37 17=1, N7 
WR I TE (61, 6700) 

6700 FORMAT (170H */*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/* 

1 /it/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/*/) 
WRITF(61 ,4)V0,V1,V2,V3,Y0,Y1,Y?,Y3 

4 F0RMAT(////10X,5HV0= F15.10»10X,5HV1= Fl 5 . 10 , 1 OX , 5HV2= Fl5.10,l 

10X,5HV3= F15.10/10X,5HY0= Fl 5 . 1 0 , 1 OX , 5HY1= Fl 5 .in , 1 OX , 5HY2= Fl 

25.10,10X,FHY3= El 6 • 10 ) 

SMU1= (Vl-VO) /Y1 
SMU?= ( V2-V1 ) /(Y2-Y1 ) 

SMU3= ( V3-V2 ) / C Y3-Y2 ) 

AK1=V1 /VO 
AK2=V2/V0 

IF( ABS(V2— VI )-l .0F-09)4112»4117,4ni 
4112 D2=9.9R9999999E206 
GO TO 4113 

4111 D2=l • 0— SMU1 /SMU2 
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At 1 3 I F ( ABS ( V3-V2 ) — 1 .0F-09J411 A. At 1 A* All 5 
AHA D3=9.999O9999QE206 
GO TO 4116 

4115 03= 1 • 0-SMU1 /SMU3 

4116 GA 1 = 1 • 0— AKl **2 
GA2=1.0-AK2**2 

RFT A= < AKl /(AK 1+1.0) )*( (Y2-Y1 )/Yl ) 

GAMA=(AK2+1.0)*( (Y2-Y1)/Y1) 

CALL SSWT CH ( 2 * J ) 

GO TO (4901 ,4902) >J 
4901 WR 1 TF ( 61 * 7 ) SMU1 »SMU2 >SMU3 

7 F0RMAT(/25X.6HMU1= FI 7 . 1 0 , 5X ♦ 6HMU2 = FI 7 . 1 0 » 5 X , 6HMU3= F17.10) 
WRITE (61 *28 1D2.AK1 *GA1 »D3 »AK2 »GA2 *GAMA ,RFTA 
28 FORMAT (2FX*6HD2= F17. 10 *5X ,6HK1= FI 7 . 10 » 5X , 6H0FLT1 =F17 . 10/2 5X » 

16HD3= F17.10*5X*6HK2= FI 7 . 1 0 » 5X , 6HDFLT 2=F 1 7 . 1 0/2 5X . 6HG AMMA =Fl 7 


4902 

8 

6005 

6006 


9 

1411 

1412 

1413 
14 

1611 

1612 

1613 

17 


173 


2.10»5X*6HRETA= E17.10/) 

I F ( ABS ( SMU1— SMU? ) — 1 • OF— 07 ) 9 » 9 » R 
I F ( APS ( SMIJ2— SMU3 ) — 1 ,0F-07)11 »11 *6005 
WRITE(61 *6006) 


FORMAT ( /47X*27H**F0CUS 
COMP= 1 2.0 
GO TO 6010 
I F ( SMU 1 ) 1 3 * 1 41 1 » 1 5 
WR I TE ( 61 * 1 41 2 ) 

FROM 

SFCOND 

LAYER**) 

FORMAT ( / 47X . 27H**F0CUS 
WR I TE ( 61 * 1 41 3 ) 

FROM 

SECOND 

layfr**) 

FORMAT ( /5?X » 1 8H**** NO 
GO TO 14 

IF(SMU3) 1611 .1611 »17 
WR I TE ( 61 » 1 61 2 ) 

FOCUS 

****) 


FORMAT ( /47X * 26H**FOCUS 
WR I TE ( 61 » 1 61 3 ) 

FROM 

THIRD 

LAYER**) 

FORMAT ( /52X.1 RH**** NO 
GO TO 16 

FOCUS 

**** ) 



WR I TE ( 61 * 1612 ) 

CN1=V0/(SMU3*Y2) 

COT GTS= SORT ( CN1 ) 

WRITE (61*173 )CN1 »COTGTS 

FORMAT (33Xt6 HQ = E 1 7 . 1 0 . 5X ,6HCTST= E17.10) 

XF=4.*Y2*COTGTS 

CN2=V0/(SMU3*Y2+V0 ) 

THET S = ACOS ( SQRT ( CN2 ) ) 

THETM= ( ACOS ( V0/V3 ) ) *57 . 2957795 1 
SCNDD=(8.*Y2*( 1.+CN1 )**2)*COTGTS 
0MEGE=(4.*C0TGTS) /(XF*SCNDD) 

THOS=THETS*57. 29577951 
A I L = ABS ( OMEGE ) *XF*XF 


WRITE (61 *2078 ) THOS »THFTM *XF » SCNOO *OMFGF ♦ A I L 
2078 FORMAT ( 33X »6HTHOS= FI 7 . 1 0 » 5X , 6HTHfTM=F1 7 . 10 /33X » 6HXF= F17.10.5X, 

16HSCNDD = F17.10/33X*6H0MEGA = F17.10»5X»6HAIL= FI 7. 10//) 

GO TO 16 

13 I F ( SMU3 )1666*1666»19 
1666 WR I TE ( 6 1 , 1 41 2 ) 

WRITE ( 61*1413 ) 

WR I TE ( 61 * 1612 ) 
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WR I TE (61*1613) 

GO TO 16 

15 I F ( SMU3— SMU1 ) 1 666 *1666 *19 

19 CN4=I (V2-V0)*CY3-Y2n/( (V3-V2>*Y2> 

CN3=2.*CN4 

CN5=CN4**2 

CN6=( V2/V0 )**2 

CN7=C1.0-CN4)**2 

CN8=< ( CN3—CN5 ) / ( CN6— CN7 ) ) 

WR I TE ( 61 *141 2 l 
WR I TE (61*1413) 

WR I TE ( 61 * 1612 ) 

Q=CN8/(1.-CN8) 

FC=1 •-D3/SQRT ( 1»+GA2*Q) 

WR I TE C 61 *241 ) FC 
CTST=SQRT (0) 

WR ITE(61»'173)Q»CTST 
C0STHS=S0RT(CN8) 

THETS*ACOS(COSTHS ) 

THETM=(ACOS(VO/V2))*57. 29577951 
CN9=( ( V2+V0 ) / ( VO— V2 ) ) 

CN10=CN9*(CN7-1.0> 

CN1 1=SQRT (CN10) 

XF=2#*Y2*CN11 

SCNDO* ( 2. 0*V0* < 1 .-AK 2**2 )* ( 1 . /SMU3-1 • /SMU1 1) / ( FOSTHS* C 1 ,-AK2**?*CN 
18 ) **1 #5 ) 

0MEGE»(4.*<C0STHS/STN(THFTS) M/(XF*SFNnr>) 

THOS=THETS*57. 29577951 
AIL=A8S(0MEGE)*XF*XF 

WR I TE (61*2050 ) THOS *THFTM . XF .SCNDD.OMFGF* AIL 
2050 FORMAT ( 33X »6HTH0S= F17. 10*5X,6HTWFTM=F17.10/33X*6HXF= F17.10»5X* 

16HSCNDD=F17.10/33X»6H0ME6A=F17.10*5X*6HA1L= E17.10//) 

GO TO 16 
11 C0MP=35.0 

IF (V3-V0) 6010. 6010*4513 
4513 THETM=( ACOS(VO/V3) ) *57 • 29577951 
6010 IF(SMU1)20. 21*22 
21 IF(SMU2)1 03 0*1030*23 
1030 WR ITE ( 61 ,1613 ) 

GO TO 10 

23 CN12=V0/(SMU2*Y1) 

FC=CN12-(Y2— Y1)/(Y1*(AK2-1. ) ) 

WRITE(61»241)FC 
C0TGTS=5QRT (CN12) 

WRITE(61*173)CN12»C0TGTS 
XF=4.*Y1*C0TGTS 
THETS=ATAN( l.O/COTGTS) 

THETM=( ACOS(VO/V2) >*57.29577951 
SCN0D=(4.*Y1*(1.+CN12)**2)/C0TGTS 
OMEGE= ( 4. *COTGTS ) / ( XF*SCNDD > 

AIL=ABS( OMEGE ) *XF*XF 
GO TO 6009 

20 IF(SMU2)1020, 1020. 2009 
1020 WRI TE ( 61 , 1613 ) 

GO TO 10 
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2009 I F ( V2-V0 )101»101»24 
101 WRITE (61*1613) 

GO TO 10 

22 I F ( SMU2-SMU1 ) 1010*1010*24 
1010 WRI TE ( 61 *1613 ) 

GO TO 10 

24 CN13=1.0-D2**2 
CN14=AK1**2-D2**2 
CN15=5QRT (CN13/CN14) 

Q=CN 15 /SORT ( l.-CN 15**2) 

CTST=SQRT (Q) 

FC=1.-(D2*SQRT(1.-CN15**2) )/SQRT(l.-< AK1 *CN 1 5 ) **2 ) 

WRITE (61*241 )FC 

241 FORMAT(/33X.17HFOCAL CONDITION= E17.10) 

WRITE(61»173)Q*CTST 

THETS=ACOS(CN15) 

THETM=( AC0S(V0/V2) ) *57. 29577951 

CN1 6= 1 .0+AK1 

CN17=1 .0— AK1 

CN18=D2** 2-1.0 

CNl 9 = (CN16/CN17) *CN1 8 

CN20=SQRT (CN19) 

XF=2.*Y1*CN20 

SCNDD=(2.*V0*( 1./SMU2— 1.0/SMU1 )*( l.-AKl**2) ) / < CNl 5* ( 1 . -AK 1**2*CN1 5 
1**2 )**1.5 ) 

COTGTS = CN 1 5/SORT ( 1. -CNl 5**2 ) 

OMEGE= ( 4.*C0T GTS ) / ( XF*SCNDD ) 

A I L = ABS ( OMEGE ) *XF*XF 
6009 THOS=THETS*57. 29577951 

WRITE(61»25) THOS »THETM»XF * SC N OD * OMEGE » A I L 
IF(COMP-12. 0)16*10*16 
18 THOS=THETS*57. 29577951 

WRITE(61*25)TH0S*THFTM»XF»SCND0.0MEGE*AIL 

25 FORMAT (33X»6HTHOS= F17 . 1 0 »5X ,6HTHFTM*F1 7 . 10/33X » 6HXF= F17.10.5X* 

16HSCNDD=F17.10/33X*6H0MEGA=E17.10*5X»6HAIL= E17.10//) 

GO TO 16 

10 IF(SMU3)1611*1611.1276 
1276 IF(V3-V0) 1611 . 1611*12 
12 WRI TE ( 6 1 , 27 ) 

27 FORMAT (//47X»26H**FOCUS FROM THIRD LAYER**) 

IF( SMU3-SMU2) 29.29.1031 
1031 IF(SMU2)38 *1100.38 
1100 D2=9.999999999E206 

38 IF(D2 )39, 43*44 

43 WRITE(61.45)SMU1 »SMU2 

45 FORMAT! /9X.29HTWO LAYER CASE WHERE MU1=MU2 E14. 5 * 1 OX * E14. 5 ) 

GO TO 16 

39 WRITE(61.6245 ) 

6245 FORMAT(57X. 11HSUB DOMAIN A) 

CN21=1.0-GA1/GA2 

CN22=D2/(D2-D3) 

S1=(1.0/CN21)*CN22 
CN23=-1.0/( D2— D3 ) 

CN24=GA1/GA2 

48 CN25=(1,0-CN24)*S1**2+CN24 
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CN26=SQRT ( CN25 ) 

S2=CN23+CN22*(S1/CN26) 

CALL SSWTCH C 1 ♦ J ) 

GO TO (6200.6201 ) , J 

6200 WR I TE (61.2191 )S1 

2191 FORMAT ( 40X»29H*ITERATION FOR EQUATION 500= E17.10) 

6201 CN27=ABSCS1-S2) 

IE(CN27-T0L1 J46.46.47 

47 S1=S2 
GO TO 48 

46 Q=(1.0/GA2>*( (1.0/S2**2)-1.0> 

FC=l.-( (02>/SQRT(l.+GAl*Q> )+( ( D2-D3 ) /SORT ( 1 .+GA?*Q J J 

WRITE(61,241 JFC 

CTST=SQRT (Q J 

WRITE(61, 173)0. CTST 

CN28=1 .O/CTST 

THET5=ATAN(CN28 ) 

THETM= ( ACOS ( V0/V3 ) >*57.29577951 
CN29=SQRT ( 1 • 0+GAl*Q ) 

CN30=SQRT(1.0+GA2*Q) 

CN31=D2-D3 

XE= ( ( 2 .*Y1 ) / ( AK1— 1 .0) ) *( 1 .O/CTST ) * ( 1 • 0— D2*CN29+CN31*CN30 ) 

SCNDD=( (?.*V0*(1.+Q)**?)/(SMU1*CTST) )*( (D2-D3 >*GA?/(1.+GA?*Q)**1.5 
1— D2*GA 1/(1 .+GA1*Q ) **1 .5 ) 

0MEGE=4.*CTST/(XF*SCN0D) 

A I L=ABS ( OMEGE ) *XF*XF 
GO TO 18 

44 IF(D2—1.0) 49.502 .51 
502 IF(GAMA-1 .0)52,53.54 
52 Sl=l • 0— GAMA* ( SMU2/SMU3 ) 

55 S2= ( ( 1.0-GAMA)*S 1**0. 3333333333 )+GAMA* Cl. 0-SMU2/SMU3) 

CALL SSWTCH ( 1 , J ) 

GO TO (6202,6203) ,J 

6202 WRITE (61,2091 )S1 

2091 FORMAT ( 40X » 29H* ITERATION FOR EQUATION 750= E17.10) 

6203 CN33=ABS( SI— S2 ) 

I F(CN33-T0L1 >56.56,57 
57 S1=S2 
GO TO 55 

56 Q=(S2**0. 6666666667-1. 0)/GA2 
CTST=SQRT (Q ) 

FC= 1 . /SMU2+ ( 1 . /SMU3-1 ./SMU2 ) /SORT ( 1 .+GA2*Q ) - ( Y1*Q ) /VO 

WR I TE ( 61,241 ) FC 

WRITE(61, 173)0, CTST 

CN34=1. O/CTST 

THETS=ATAN(CN34) 

THETM=(ACOS(VO/V3) >*57.29577951 
CN35=SQRT ( 1 . 0+GA2*Q ) 

XF=4.*YI *(!.+ ( (AK2+1 .)/?.)*( (Y2-Y1 )/Yl )* ( 1 .-SMIJ2/SMU3 ) /CN35 ) *CTST 
SCN0D=(7.*V0*(1 .+0)**?/CTST)*(GA?*(l ./SMU3-1 ./SMU2)/( l.+GA2*0)**l. 
15+2 •* Y1 /VO ) 

0MEGE=4.*CTST/(XF*SCN0D) 

A I L = ABS ( OMEGE ) *XF*XF 
GO TO 18 

54 Sl=l./( (l.-SMU2/SMU3)**3) 
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58 S2= (l.+( GAMA-1. ) *51**0. 6666666667) /(GAMA* (1.-SMU2/SMU3) > 

CALL SSWTCH ( 1 » J ) 

GO TO (6204,6205) ,J 

6204 WRITE(61,59)S1 

59 FORMAT (40X.29H* HERAT I ON FOR EQUATION 775= E17.10) 

6205 CN36=ABS(S1-S2) 

I F ( CN36-TOL 1)60,60,61 

61 S1=S2 

GO TO 58 

60 Q=( l./GA2)*( ( l./S2**0. 6666666667 )-l. ) 

FC=1 » — ( (D2)/SQRT(1.+GA1*Q) )+( ( D2-D3 ) /SORT ( 1 ,+GA2*Q ) ) 

WRITE ( 61,241 )FC 
CTST=SQRT (Q ) 

WRITE (61,173)0»CTST 
CN37* 1 • /CTST 
THET S=ATAN ( CN37 ) 

THETM=( AC0S(V0/V3) )*57. 29577951 
CN38=SQRT(1.+GA2*Q) 

XF=4.*Y1*(1 .+( ( AK2+1. )/2. )*( (Y2-Y1 ) /Y1 ) * ( 1 . -SMU2/SMU3 ) /CN38 ) *CTST 
SCNDD«( 2.*V0*( l.+Q)**2/CTST )*(GA2*( 1./SMU3-1./SMU? )/( l.+GA2*0)**l. 
15+2.*Y1/V0> 

0MEGE=4.*CTST/(XF*SCNDD> 

AIL=ABS(OMEGE)*XF*XF 
GO TO 18 

53 Q=( 1.-SMU2/SMU3 ) **0.6666666667-1.0 
CTST=SQRT (Q ) 

WRITE(61,173)Q,CTST 
CN39= 1 . /CTST 
THETS=ATAN( CN39 ) 

THETM= ( ACOS ( V0/V3 ) )*57. 29577951 
CN40= SORT ( 1 .+GA2*Q ) 

XF=4.*Y1*(1 .+( ( AK2+1. )/2. )*( ( Y2-Y1 ) /Y1 ) * < 1 .-SMU2/SMU3 ) /CN40 ) *CTST 
SCNDD=(2.*V0*(1.+Q)**2/CTST)*(GA2*(1./SMU3-1./SMU?)/(1.+GA2*Q)**1. 
15+2.*Y1/V0) 

0MEGE=4.*CTST/ ( XF*SCNDD ) 

AIL=ABS(OMEGE)*XF*XF 
GO TO 18 

49 WR I TE ( 61 ,6246 ) 

6246 FORMAT(57X»l 1HSUBDOMAI N B) 

Sl=l.+( (D3-1.)/(1.-D2*(1.-GA1/GA2) ) ) 

62 CN41= ( l.-GAl/GA2)+( GA1 /GA2 ) *S1**2 
CN42=SQRT(CN41 ) 

S2=D3-D2+D2*S1/CN42 

CALL SSWTCH (1,J) 

GO TO (6206.6207) ,J 

6206 WR I TE ( 61 ,63 ) SI 

63 FORMAT ( 40X , 29H* I TERATION FOR EQUATION 600= E17.10) 

6207 CN43=ABS(S1-S2) 

IF(CN43-T0L1 )64,64,65 

65 SI =S2 
GO TO 62 

64 Q=(S2**2-1. )/GA2 

FC= 1 • — ( (D2)/SQRT(1.+GA1*Q) )+( ( 02-D3 ) /SORT ( 1 ,+GA2*Q ) ) 

WRITE(61,241 )FC 
CTST=SQRT(Q) 
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WRITEI61 ,173)Q»CTST 
CN44=1./CTST 
THETS=ATAN( CN44 ) 

THETM=(AC0S(V0/V3> >*57.29577951 

XF=(?.*Y1/CAK1-1. ) )*( 1 • /CTST ) * ( 1 •— D2*SQRT ( 1 .+GA1*Q ) + ( D2-D^ ) *SQRT ( 1 
l.+GA2*Q) ) 

SCNDD= ( (2.*V0*< l.+Q)**2)/(SMUl*CTST) )*< ( D2~D3 ) *GA2/ < 1 • +GA 2*Q > **1 • 5 
l-D2*GAl/( l.+GAl*Q)**1.5) 

0MEGE=4.*CTST/(XF*SCNDD) 

A I L=ABS(OMEGF)*XF*XF 
GO TO 18 

51 IF(D2-Y2/Yl >66.67.68 

67 CN45=(l.-( ( Y2— Y1 ) /Y2 ) * ( SMU2/SMU3 ) )**2 
Q=( 1./GA1 ) * ( 1 • /CN45— 1 ■ ) 

FC=1 •— ( D2/SORT ( 1 • +GA 1*0) ) + ( ( D2-D3 ) /SQRT ( 1 .+GA2*Q ) > 

WRITE(61*241)FC 
CTST=SQRT (0 ) 

WRITE(61»173)G»CTST 
CN46=1 • /CTST 
THETS=ATAN ( CN46 ) 

THETM= ( ACOS ( V0/V3 > >*57.29577951 
XF=2.*( AK1+1. )*Y2*SQRT(Q/(1.+GA1*Q> ) 

SCNDD=-(2.*V0*D2*GA1*(1.+Q)**2)/(SMU1*CTST*( l.+GAl*Q>**1.5> 
0MEGE=4.*CTST/(XF*SCNDD) 

AI L= ABS ( OMEGE J *XF*XF 
GO TO 18 

66 WRITE(61.6247) 

6247 FORMAT ( 57X * 11HSUBD0MAIN C) 

Sl=1.0 

CN47* ( 1 GA1/GA2 ) 

69 CN48= ( ( 02/ ( S1+D2-D3 ) )**2 J-GA1/GA2 
S2= SORT ( CN47 /CN48 ) 

CALL SSWTCH ( 1 » J ) 

GO TO (6208.6209) »J 

6208 WR I TE (61.70)51 

70 FORMAT ( 40X »29H* ITERAT ION FOR EQUATION 800= E17.10) 

6209 CN49=ABS ( SI— S2 ) 

IF( CN49— T0L1 >71.71.7 2 
72 SI =S2 
GO TO 69 

71 Q=(S2 **2—1.0) /GA2 

FC=1 •— ( ( 02 ) /SORT ( 1 .+GA1*Q) )+( ( 02-03 ) /SORT ( 1 .+GA2*0 > ) 

WR I TE ( 61 *241 ) FC 
CTST=SQRT (Q ) 

WRITE(61, 173)0. CTST 
CN50=1 • /CTST 
THETS=ATAN(CN50 ) 

THETM=(AC0S(V0/V3) >*57.29577951 

XF=(2.*Y1/(AK1-1. ) )*(1./CTST)*(1.-02*SQRT(1.+GA1*0)+(02— D3)*SQRT(1 
l.+GA2*Q) ) 

SCNDD= ( (?.*V0*(l.+0)**2) /(SMU1*CTST) )*( ( D2-03 ) *GA2/ ( 1 .+GA2*0 >**1 . 5 
1— D2*GAl/( l.+GAl*0)**l,5 ) 

0MEGE=4.*CTST/( XF*SCNDD ) 

AIL=ABS(OMEGE)*XF*XF 
GO TO 18 


25 



68 IF(SMU2)73»7431 *73 
73 WR I TE ( 61 *6248 ) 

6248 FORMAT ( 57X»1 1HSUBDOMAI N D) 

Sl=< 1 • / < 1.-GA2/GA1 ) )*( (D2-D3J/D2) 

75 CN1 0= ( 1.-GA2/GA1 )*S1**2+GA2/GA1 

S2 = l./D2 + ( (D2-D3) /D2 >*< S1/SQRT(CN10) ) 

CALL SSWTCH ( 1 * J ) 

GO TO (6210*6211 ) *J 

6210 WR I TE ( 6 1 *147 ) SI 

147 FORMAT(40X»30H* ITERATION FOR EQUATION 1000= E17.10) 

6211 IF(ABS(S1-S2)-T0L1 >77*77*78 

78 S1=S2 

GO TO 75 

77 Q=(1./GA1)*(1./S2**2-1.) 

FC=1 • - ( ( D2 ) /SORT (l.+GA 1*0) )+( < D2-D3 ) /SORT ( 1 .+GA2*Q ) ) 

WR I TE ( 61 ,241 ) FC 
CTST=SQRT (Q ) 

WRITE(61,173)Q»CTST 
CN1 1= 1 • /CTST 
THET S = ATAN(CN11 ) 

THETM=(AC0S(V0/V3) >*57.29577951 

XF=(2.*Y1/(AK1-1 ) )*< l./CTST)*< l.-D2*SQRT(l.+GAl*Q)+(D2-D3)*SQRT( 1. 
1+GA2*Q ) ) 

SCNDD=( ( 2 . *V0* ( 1 « +Q ) **2 ) / (SMU1*CTST) ) *( ( D2-D3 > *GA2 / ( 1 • +GA2*Q ) **1 . 5 
1— D2*GA 1/(1 »+GA 1*Q ) ** 1 • 5 ) 

0MEGE=4.*CTST/( XF*SCNDD ) 

AIL=ABS(OMEGE)*XF*XF 
GO TO 18 

7431 IF(D3-BETA)79, 80*81 

79 S1=(2./(3.*BETA) >+( ( BET A-D3 > /BETA > **1 . 5 

76 S2=( l. + (BETA-D3)*Sl **0.3333333333 ) /BET A 
CALL SSWTCH (1,J) 

GO TO (6212,6213) »J 

6212 WR ITE(61,92)S1 

92 FORMAT ( 40 X»30H* ITERATION FOR EQUATION 1175= E17.10) 

6213 IF(ABS( S1-S2 )-TOLl >145,145,146 
146 SI =S2 

GO TO 76 

145 Q= ( 1 * /GA1 )*( l./S2**0. 6666666666-1. > 

FC=1.-D3/SQRT( 1 .+GA1*Q)+(Q*GA1*AK1*( Y2-Y1 ) ) / ( Yl* ( AK1 + 1 . )* ( l.+GAl*Q 
1 > ** 1 • 5 > 

WR ITE(61,241)FC 
CTST=SQRT (Q) 

WRITE(61,173)Q»CTST 
CN1 2= 1 • /CTST 
THETS=ATAN(CN12 ) 

THETM=(AC0S(V0/V3) >*57.29577951 

XF=(2.*Y1/CTST)*( ( 1 • -D 3* SORT ( 1 . +GA 1*Q ) )/(AKl-l. >+( (Y2-Y1 ) /Yl )*(AK1 
1*Q/SQRT ( 1 .+GA1*Q ) ) > 

SCNOD=( 2. *V0* ( 1 • +Q ) **2/ (SMU1*CTST) )*( (SMU1*( Y2-Y1 ) * ( 2. *AK 1**2-AK 1* 
1*2*GA1*Q ) )/(Vl*(l.+GAl*Q)**2.5 )-(D3*GAl/ (l.+GAl*Q)**1.5 ) ) 
0MEGE=4.*CTST/(XF*SCNDD> 

AIL=ABS(OMEGE)*XF*XF 
GO TO 18 

80 Q= ( BET A**0. 6666666667-1. 0 ) /GA1 
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CTST=SQRT(Q) 

WRITE(6l»173)Q,CTST 

CN13=1#/CTST 

THETS=ATAN(CN13) 

THETM=( ACOS (V0/V3) >*57.29577951 

XF=(2.*Y1/CTST)*( (l.-D3*SQRT(l.+GAl*Q) ) / ( AK 1-1 . 0 ) + ( ( Y2-Y1 ) / Y1 ) * ( AK 
11*0 /SORT ( l.+GAl*Q) ) ) 

SCNDD=(6.*Yl*(AKl+l. )*( l.+Q>**2) /(SQRT(Q)*( l.+GAl*Q) ) 

OMEGE= ( 4.*CTST ) / ( XF*SCNDD ) 

AIL=ABS(OMEGE)*XF*XF 
GO TO 18 

81 S1=(D3-BFTA)**3+3.*BETA 

82 S2=BETA+(D3-BETA)*S1**0. 6666666667 
CALL SSWTCHCl.J) 

GO TO (6214,6215) ,J 

6214 WR I TE ( 61 » 83 ) SI 

83 FORMAT ( 40X ♦ 30H* ITERATION FOR FOUATION 1125= F17.10) 

6215 1 F ( ABS ( SI— S2 )— T0L1 184.84*85 

85 S1=S2 

GO TO 82 

84 Q=(S2**0. 66666666667-1. 0)/GAl 

FC=1 •— D3/SQRT ( 1 .+GA1*Q )+(Q*GAl*AKl*(Y2-Yl ) > / ( Yl* ( AK1+1 . )* ( 1 .+GA1*Q 
1 ) **1 • 5 ) 

WRITE (61,241 )FC 
CTST=SORT (O ) 

WRITE(61, 173)0, CTST 
CN14=1 .O/CTST 
THETS=ATAN(CN14) 

THETM= ( ACOS ( V0/V3 ) >*57.29577951 

XF= ( 2.*Y1 /CTST ) * ( (l.-D3*SQRT( l.+GAl*Q) ) / ( AK1-1. )+( ( Y2-Y1 ) /Y1)*(AK1 
1*0 /SORT ( 1 .+GA1*0 ) ) ) 

SCND0=(2.*V0*(l.+0)**2/(SMUl*CTST) )*( ( SMU1* ( Y2-Y1 ) * ( 2. *AK 1**2-A<1* 
1*2*GA1*0) )/(Vl*(l.+GAl*0)**2.5)-(D3*GAl / ( 1 .+GA1 *0 ) **1 . 5 ) ) 
0MEGE=4.*CTST/(XF*SCNDD) 

A I L= ABS ( OMEGE ) *XF*XF 
GO TO 18 

29 I F ( D2 ) 87 ,87 ,86 
87 WR I TE ( 61 ,6249 ) 

6249 FORMAT ( 57X , 1 1HSUBD0MA I N A) 

WR I TE (61, 1613) 

GO TO 16 

86 IF(D2— 1.0)89,90,91 

90 CSTR= ( ( l.-GAMA ) /3. ) **3— ( GAMA/2,* ( SMU2/SMU3— 1 . ) )**2 
WR I TE (61,9088) CSTR 
9088 FORMAT ( 33X »6HCSTAR = E17 • 10 ) 

IF(CSTR)87,93 ,94 
93 Q=( ( (l.-GAMA)/3.)-l.)*l./GA2 
CTST=SQRT (Q ) 

WRITE(61, 173)0, CTST 
CN15=1*/CTST 

XF=4.*Y1*CTST*( l.+( AK2+1. ) /2.*( Y2-Y1 ) /Yl*( 1 .-SMU2/SMU3 ) /SORT ( 1 .+GA 
12 * 0 ) ) 

THETS=ATAN(CN15 ) 

THETM= ( ACOS ( V0/V3 ) >*57.29577951 

SCNOD= ( 2.*V0* ( 1 ,+Q )**2/CTST)*(GA2*( 1./SMU3— 1./SMU2 ) /( 1 ,+G A2*Q ) **1 • 
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15+2.*Y1/V0> 

0MEGE=4.*CTST/(XF*SCNDD) 

AIL=ABS(OMEGE)*XF*XF 
GO TO 18 

94 J2P=-1 

SI =( GAMA/ (GAMA-1.0) ) * ( SMU2/SMU3-1 . 0 ) 

95 S2=(S1**3+GAMA*(SMU2/SMU3-1.0) )/(1.0-GAMA) 

IF( ABS(S1-S2 )-TOLl ) 96,96. 97 

97 SI =S2 

CALL SSWT CH ( 1 , J ) 

GO TO (6216,6217 ) ,J 

6216 WR I TE ( 61 ,98 ) SI 

98 FORMAT (40X»30H*ITERATI ON FOR EQUATION 1375= E17.10) 

6217 GO TO 95 

96 Q=(S2**2-1.0)/GA2 

FC=l.-( ( 02) /SORT ( l.+GAl*Q) )+( ( D2-D3 ) /SORT ( 1 ,+GA2*Q ) ) 

WR ITE(61»241)FC 
CTST=SORT ( Q ) 

WRITE(61,173)Q,CTST 
CN25= 1 • /CTST 
THETS=ATAN(CN25 ) 

THETM= ( AC OS ( V0/V3 ) ) *57 . 2957795 1 

XF=4.*Y1*(1 .+( ( AK2+1 • ) / 2 • ) * ( ( Y2-Y1 ) /Y1 ) * ( 1 . -SMU2 /SMU3 ) /SORT ( 1.+GA2 
1*0) ) * ( CTS T ) 

SCNDD=(2.*V0*(1 ,+Q)**2/CTST)*(GA2*(l«/SMU3-l./SMU2 )/( 1 .+GA 2*Q > **1 . 
15+2.*Y1/V0) 

0MEGE=4.*CTST/(XF*SCNDD) 

AIL=ABS(OMEGF)*XF*XF 
J1P = 0 

100 TH0S=THETS*57. 29577951 

WRITE(61,25) THOS ,THETM,XF ,SCNDD. OMEGE , A I L 
I F ( J1P ) 16 » 1 02 » 16 

102 I F ( J2P) 103,104,105 

103 J1P=4 

Sl=( l.-GAMA)*(SMU2/SMU3) 

106 S2=( 1.0-GAMA)*S1**0.3333333333+GAMA*( 1 .0-SMU2/SMU3) 

CALL SSWT CH ( 1 , J ) 

GO TO (6218,6219) ,J 

6218 WR I TE ( 61 , 1 07 ) S 1 

107 FORMAT (40X, 3 OH* ITERATION FOR EQUATION 1376= E17.10) 

6219 IF(ABS(S1-S2)— TOL1 ) 108, 108, 109 
109 S1=S2 

GO TO 106 

108 Q=(S2**0. 6666666667-1. 0)/GA2 

FC= 1 • - ( (D2 ) /SORT ( l.+GAl*Q) )+( ( D2-D3 ) /SQRT ( 1 .+GA2*Q ) ) 

WRITEI61.241 )FC 
CTST=SQRT (Q ) 

WR ITE(61»173)Q,CTST 
CN3 1=1. /CTST 
THETS=ATAN(CN31 ) 

THETM=(ACOS(VO/V3) ) *57 . 2957795 1 

XF=4.*Y1* ( 1 .+ ( ( AK2+1. ) / 2 - ) * ( < Y2-Y1 ) /Y1 ) * ( 1 . -SMU2/SMU3 ) /SQR T ( 1.+GA2 
1*Q ) ) * ( CTST ) 

SCNDD=(2.*V0*( l.+Q)**2/CTST)*(GA2*(l./SMU3-l./SMU2 )/( l.+GA2*Q)**l. 
15+2.*Y1/V0) 


28 


0MFGE=4.*CTST/(XF*SCNDD) 

AIL=ABS(OMEGF)*XF*XF 
GO TO 100 
89 WRI TE {61*6251 ) 

6251 FORMAT ( 57X » 1 1HSUBDOMA IN B) 

CC01= ( 02**2 /( 1. 0-GA1 /GA2 )) **0.3333333333 
CC02 = ( (D2-D3)**2/(GA2/GA1-1.0> ) **0 . 3333333333 
CE=CC01-CC02 
WRITE(61,8991)CE 
8991 FORMAT (/33X,6HCE= E17.10) 

IFCCE-1. 0)110. 111*112 

110 WR I TE (61*113) 

113 FORMAT ( /52X » 18H**** NO FOCUS **** ) 

GO TO 16 

111 Q=( (D2*( 1.-GA1/GA2) ) **0 . 6666666667-1 . 0 ) /GA1 
CTST=SQRT (0) 

WRITF(61, 173)0, CTST 
CN44=1./CTST 
THETS=ATAN ( CN44 ) 

THETM=(ACOS(VO/V3) 1*57.29577951 

XF=(2.*Y1)/(AK1-1.0)*1./CTST*(1.-D2*SQRT(1.+GA1*Q)+(D2-D3)*SQRT(1. 
1+GA2*Q ) ) 

SCNDD=( (2.*V0*(1.+Q)**2)/(SMU1*CTST) ) * ( ( D2-D3 > *GA2/ ( 1 .+GA2*Q > **1 . 5 
1— D2*GA 1/ { 1 .+GA1 *0 ) **1.5 ) 

0MEGE=4.*CTST/ ( XF*SCNDD ) 

AIL=ABS(OMEGE)*XF*XF 
GO TO 18 

112 J1P = 0 
J2P=0 

SI = ( D2 /SORT ( 1 •— GA1 /GA2 ) — 1 • ) / ( D2—D3 ) 

114 CN45=SQRT ( <1.-GA1/GA2)*S1**2+GA1/GA2) 

S2= ( D2/ ( 02— D 3 ) )*S1/CN45— l./( D2— D3 ) 

CALL SSWTCH ( 1 * J ) 

GO TO (6220,6221) ,J 

6220 WRITE(61,115)S1 

115 FORMAT ( 40X » 3 OH* 1 TERATI ON FOR EQUATION 1275= E17.10) 

6221 IF(ABS(S1-S2)-T0L1 ) 116, 116, 117 
117 S1=S2 

GO TO 114 

116 Q=(1./GA2)*(1./S2**2-1.0) 

FC=l.-( (02) /SQRT(1.+GA1*Q) )+( (02-03 ) /SORT ( 1 ,+GA2*0 )) 

WR I TE ( 61 ,241 ) FC 
CTST=SQRT (Q ) 

WRITE(61, 173)0, CTST 
CN36=1./CTST 
THETS=ATAN ( CN36 ) 

THETM=(ACOS(VO/V3) )*57. 29577951 

XF=(2.*Y1/(AK1-1.) )*1. /CTST* (l.-02*SQRT(l. +GAl*0)+( 02-03 )*SQRT(1.+ 
1GA2*0 ) ) 

SCNDD= ( ( 2 .*V0* ( 1 .+0 )**2 ) / < SMU1*CTST ) ) * ( ( 02-03 )*GA?/ ( ) .+GA2*0 >**1 .5 
1— D2*GA 1/(1 .+GA 1*0 ) **1 • 5 ) 

0MEGE=4.*CT ST/ ( XF*SCNDD ) 

AI L=ABS( OMEGE )*XF*XF 
GO TO 100 
104 J1P=4 
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S1=1.-(1.-D3)/(1.-(D2-03)*(GA2/GA1-1. ) ) 

118 CN55=SQRT ( ( 1 .-GA2/GA1 )+GA2/GAl*Sl**2 ) 

S2=D2-( (D2-D3)*S1 )/CN55 

CALL SSWT CH ( 1 * J ) 

GO TO (6222,6223 ) ,J 

6222 WR I TE ( 6 1 , 1 19 ) SI 

119 FORMAT (40X.30H* ITERATION FOR EQUATION 1276= F17.10) 

6223 IF (ABS( SI -S2 )-TOLl ) 120, 120,121 
121 S1 = S2 

GO TO 118 

120 Q=(S2**2-1.)/GA1 

FC=l.-( (D2)/SQRT ( l.+GAl*Q> ) + ( ( D2-D3 ) /SORT ( 1 ,+GA2*Q ) ) 

WRITE(61 ,241 ) FC 
CTST=S0RT (Q) 

WRI TE ( 61 , 173 ) Q , CTST 
CN56=1 »/CTST 
THETS*ATAN(CN56) 

THETM=(ACOS( V0/V3) ) *57 . 2957795 1 

XF=(2.*Y1/(AK1-1. ) )*1./CTST*( l.-D2*SQRT( l.+GAl*Q)+(D2-D3)*SQRT(l.+ 
1GA2*Q ) ) 

SCNDD= ( < 2,*V0*( 1 ,+Q ) **2 ) / (SMU1*CTST ) )* ( (02-03 >*GA2/( 1 . +GA2*Q ) **1 . 5 
1— D2*GA1/(1 ,+GA 1*Q ) ** 1 • 5 ) 

0MEGE=4.*CTST/(XF*SCNDD) 

A I L = ABS ( OMEGE ) *X F*XF 
GO TO 100 
91 TST1=Y2/Y1 

WRITE(61,9111 )TST1 
9111 FORMAT (/33X,6HY2/Y1=E17,10) 

I F ( D2-TST1 >122,123,124 

124 WR I TE ( 61 , 6253 ) 

6253 FORMAT ( 57X , 1 1 HSUBDOMA I N D) 

I F ( SMU2 ) 125,126,125 

126 WR I TE ( 61 » 1 27 ) 

127 FORMAT! /52X,18H**** NO FOCUS ****) 

GO TO 16 

125 S1=D2-SQRT (GA1/GA2)*(D2-D3) 

128 CN65=SQRT ( ( 1 .-GA2/GA1 )+GA2/GAl*Sl**2 ) 

S2=D2-((D2-D3)*S1)/CN65 

CALL SSWTCH ( 1 , J ) 

GO TO (6224,6225 ) ,J 

6224 WR I TE ( 61 , 1 29 ) SI 

129 FORMAT (40X ,3 OH* I TERATION FOR EQUATION 1600= F17.10) 

6225 I F( ABS (SI -S2 )-TOLl >130,130,131 
131 SI =S2 

GO TO 128 

130 Q=(S2** 2-1,0) /GA1 

FC=l,-( (02 )/SQRT ( 1 »+GAl*Q ) ) + ( ( 02-03 ) /SORT ( 1 ,+GA2*Q ) ) 

WRITE (61,241 )FC 
CTST=SQRT (Q) 

WRITE(61,173)Q,CTST 
CN87= 1 , /CTST 
THETS=ATAN(CN87> 

THETM= ( A COS ( V0/V3 ) >*57.29577951 

XF=(2.*Y1/(AK1-1. ) )* 1, /CTST* ( 1. -02*SQRT(1.+GA1*Q)+ (02-03 )*SQRT(1.+ 
1GA2*Q ) ) 
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SCNDD=( (2.*V0*( l.+Q)**2)/(SMUl*CTST) >*( ( D2-D3 ) *GA2/ < 1 .+GA2*Q ) **1 . 5 
1-D2*GA1/(1.+GA1*0)**1.5) 

0MEGE=4.*CTST/(XF*SCNDD) 

A I L=ABS ( OMEGE ) *XF*XF 
GO TO 18 

123 COl=SMU2/SMU3 
C02 = Y2 / ( Y3-Y2 ) 

I F(C01-C02) 132, 132,133 

133 WRITE ( 61 » 127 ) 

GO TO 16 

132 C03=(l.-( (Y2-Y1 >/Y2)*SMU2/SMU3)**2 
Q=1./GA1*(1./C03-1.0> 

XF=2.*(AK1+1. )*Y2*SQRT(Q/(1.+GA1*Q) ) 

CTST=SQRT (0 ) 

FC=l.-( ( D2 ) /SORT ( 1 .+GA1*Q ) )+( ( D2-D3 ) /SORT ( 1 »+G A2*Q ) ) 

WRI TE( 61 .241 ) FC 
WRITE(61,173)Q»CTST 
CN95=1.0/CTST 
THETS=ATAN(CN95 ) 

THETM=(AC0S(V0/V3) >*57.29577951 

SCNDO*- ( 2 .*V0*D2*GA1* ( 1 ,+Q )**2 ) / ( SMU1*CTST* ( 1 .+GA1*Q ) **1 . 5 ) 
0MEGE=4.*CTST/(XF*SCN0D) 

AIL=ABS(OMEGE)*XF*XF 
GO TO 18 

122 WRI TE ( 61 *6252 ) 

6252 FORMAT ( 57X»11HSUBD0MAIN C) 

CC03= ( D2**2/ ( 1 .0— GA1 /GA2 ) >**0.3333333333 
CC04= ( ( D2— D3 ) **2/ ( 1 .0— GA2/GA1 ) >**0.3333333333 
CE= CC03+CC04 
WRITE(61,1221 )CE 
1221 FORMAT ( /33X »oHCE= E17.10) 

IF(CE-1.0)136»135»134 

134 WR I TE ( 61 * 1 27 > 

GO TO 16 

135 C04=(D2*(1.-GA1/GA2) >**0.6666666667-1.0 
Q=C04/GA1 

CTST=SQRT (Q ) 

WRITE(61»173)0»CTST 

CN75=1.0/CTST 

THETS=ATAN(CN75> 

THETM= ( ACOS ( V0/V3 ) >*57.29577951 

XF=(2.*Y1/(AK1— 1. > )*1./CTST*( l.-D2*SQRT( 1 .+GA1*Q >+ { D2-D3 > *SQRT ( 1 .+ 
1GA2*0 ) ) 

SCNDD=( (2.*V0*(1.+Q>**2) /(SMU1*CT5T) )*( ( D2-D3 > *GA?/ ( 1 .+GA2*Q > **1 .5 
1-D2*GA1/(1.+GA1*Q)**1.5) 

0MEGE=4.*CTST/(XF*SCNDD) 

A I L = ABS ( OMEGE > *XF*XF 
GO TO 18 

136 J1P=0 
J2P = 1 

S1=(D2/S0RT( l.-GAl/GA2)-1.0> /(D2-D3> 

137 C 06= SORT ( ( 1 ,-GAl /GA2 ) *5 1**2+GA 1 /GA2 ) 

S2=(D2/(D2-D3) > *S1 /C06-1 . / ( D2-D3 > 

CALL SSWT CH ( 1 » J ) 

GO TO (6226,6227) ,J 


31 



6226 WRITE C 61*138) SI 

138 FORMAT ( 40 X.30H* ITERATION FOR EQUATION 1475= E17.10) 

6227 IF( ABS( SI -S2)-T0L1 >139.139*140 

140 S1=S2 

GO TO 137 

139 Q=1./GA2*(1./S2**2-1. ) 

FC= 1 • — ( ( D2 ) /SQRT(1.+GA1*Q) )+( ( D2-D3) /SORT ( 1 .+GA2*Q ) ) 

WRITE ( 61*241 )FC 
CTST = SORT ( Q ) 

WRITE(61,173)0»CTST 
CN77= 1 • 0/CT ST 
THETS=ATAN(CN77) 

THETM=(AC0S(V0/V3) 1*57.29577951 

XF=(2.*Y1/(AK1— 1. >)*1./CTST*<1.-D2*SQRT(1.+GA1*QJ+<D2-D3»*SQRTI1.+ 
1GA2*Q> ) 

SCNDD= ( (2.*V0*(1.+Q)**2)/(SMU1*CTST) )*( ( D2-D3 ) *GA2/ ( 1 . +GA2*Q ) **1 . 5 
1-02*GA 1/ ( 1 • +GA1*Q ) **1 • 5 ) 

0MEGE=4.*CTST/(XF*SCNDD) 

AIL=ABS(0MEGE)*XF*XF 
GO TO 100 
105 J1P=4 

S1=1.+(D3-1. )/( l.+(D2-D3)*( 1.-GA2/GA1 > ) 

141 C07=SQRT( Cl.-GA2/GA1)+(GA?/GA1)*S1**2> 

S2=D2- ( ID2-D3)*S1 )/C07 

CALL SSWTCH ( 1 ♦ J ) 

GO TO (6228 .6229) »J 

6228 WR I TE ( 61 * 1 42 ) SI 

142 FORMAT ( 40X » 30H* ITERAT I ON FOR EQUATION 1476= E17.10) 

6229 I F C ABS(S1-S2 )-TOLl ) 143. 143. 144 
144 S1=S2 

GO TO 141 

143 Q=(S2**2-1. )/GAl 

FC=l.-( (D2 ) /SORT ( l.+GAl*Q) )+( ( D2-D3 ) /SORT ( 1 .+GA2*0 ) ) 

WR I TE ( 61*241 ) FC 
CTST=SQRT (Q ) 

WRITE(61»173)Q»CTST 
CN91=1 .0/CTST 
THETS= ATAN ( CN91 ) 

THETM=(ACOS( V0/V3) >*57.29577951 

XF= (2.*Y1/ ( AK1-1.0) )*1./CTST*( l.-D2*SQRT(l.+GAl*Q)+(D2-D3)*SQRT(l. 
1+GA2*Q ) > 

SCNDD= ( (2,*V0*(l.+O)**2)/(SMUl*CTST) )*( ( 02-03 > *GA2/ (1 .+GA2*Q ) **1 • 5 
1-D2*GA1/(1.+GA1*Q)**1.5 ) 

0MEGE=4.*CTST/(XF*SCNDD) 

AIL=ABS(OMEGE)*XF*XF 
GO TO 100 
16 Y3 = Y3 +DY3 
37 CONTINUE 
Y2=Y2+DY2 
36 CONTINUE 
Y1 = Y1+DY1 
35 CONTINUE 
Y0=Y0+DY0 
34 CONTINUE 
V3=V3+DV3 
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33 CONTINUE 
V2= V2+DV2 
32 CONTINUE 
V1=V1+DV1 
31 CONTINUE 
VO=VO+DVO 
30 CONTINUE 

IF( LAST) 1000, 1001, 1000 
1001 STOP 
END 

FINIS 
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APPENDIX II 


AVERAGE INTENSITY LEVEL IN THE NEIGHBORHOOD OF A FOCUS 


A 


From Reference 1, pp. 34-36, it is seen that the mean density, p^, 
over a small stretch of the x-axis is 


Rn = 




(A-l) 


where p r 2 is related to the energy output of the sound source and a) is 
a small quantity dependent on the atmospheric parameters. The quantity 
p r 2 (from Ref. 4) is equal to 


JL JL 

2* V Q ’ 


where P is the total acoustic power of the source and V Q is the propagation 
velocity of sound at the sound source. Upon choosing an initial reference 
intensity. 


T , 1n -i 2 watts 


it follows that the average intensity level in the € neighborhood of a 
focus may be given by 


I.L. =0-18 + 10 log 1Q |w| - 10 log 1G € db re 10" 12 watts/m 2 , (A-2) 


where D is the total power level of the sound source in db re 10~ 13 watts. 

Equation (A-2), then, enables one to compute the average intensity 
level quite readily. Also, it leaves one free to choose not only the 
power output of the source but also the neighborhood about the focus 
over which he wishes to know the average intensity level. It must be 
emphasized here that £ must not be chosen too large, since the express- 
ion for U) is derived for an infinitesimal neighborhood. 
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